#include "numuCCAnalysis.hxx"
#include "FiducialVolumeDefinition.hxx"
#include "Parameters.hxx"
#include "UseGlobalAltMomCorrection.hxx"
#include "StandardSystematicsNumuCC.hxx"
#include "numuCCSelection.hxx"
#include "fgd2CCSelection.hxx"
#include "CategoriesUtils.hxx"
#include "BasicUtils.hxx"

//********************************************************************
numuCCAnalysis::numuCCAnalysis(AnalysisAlgorithm* ana) : baseAnalysis(ana) {
//********************************************************************

  // Minimum accum level to save event into the output tree
  SetMinAccumCutLevelToSave(ND::params().GetParameterI("numuCCAnalysis.MinAccumLevelToSave"));
  SetMinAccumCutLevelToSaveRooTracker(ND::params().GetParameterI("numuCCAnalysis.MinAccumLevelToSaveRooTracker"));

  // Add the package version
  ND::versioning().AddPackage("numuCCAnalysis", anaUtils::GetSoftwareVersionFromPath((std::string)getenv("NUMUCCANALYSISROOT")));

  // which analysis: FGD1, FGD2 or FGDs
  _whichFGD = ND::params().GetParameterI("numuCCAnalysis.Analysis.Sample.whichFGD");

  // For FGD2, same categories with prefix "fgd2", i,e, "fgd2reaction" etc.
  ND::categ().AddStandardCategories("fgd2");
}

//********************************************************************
void numuCCAnalysis::DefineSelections(){
//********************************************************************

  // ----- Inclusive CC -----------
  if (_whichFGD==1) // FGD1
    sel().AddSelection("kTrackerNumuCC", "inclusive numuCC selection",     new numuCCSelection(false)); // true/false for forcing break
  else if (_whichFGD==2) // FGD2
    sel().AddSelection("kTrackerFGD2CC", "inclusive numuCCFGD2 selection", new fgd2CCSelection(false)); // true/false for forcing break
//  else if (_whichFGD==3) // both FGDs, in to branches
//    sel().AddSelection("kTrackerFGDsCC", "inclusive numuCCFGDs selection", new numuCCFGDsSelection(false)); // true/false for forcing break
}

//********************************************************************
void numuCCAnalysis::DefineCorrections(){
//********************************************************************

  // Some corrections are defined in baseAnalysis
  baseAnalysis::DefineCorrections();

  //------------ Corrections -----------

  // Change the main fit momentum by the muon alternate momentum, but only for P5
#if !VERSION_HAS_EQUIVALENT_MAIN_AND_ALT_FITS
  if (ND::params().GetParameterI("numuCCAnalysis.Corrections.EnableUseMuonAltMom") == 1) {
    corr().AddCorrection("altmom_corr", new UseGlobalAltMomCorrection(UseGlobalAltMomCorrection::kMuon));
  }
#endif

}

//********************************************************************
void numuCCAnalysis::DefineSystematics(){
//********************************************************************

  // Some systematics are defined in baseAnalysis
  baseAnalysis::DefineSystematics();

  if (ND::params().GetParameterI("numuCCAnalysis.Systematics.EnableMomResOptimization"))
    syst().ReplaceVariationSystematic(baseAnalysis::kMomResol, new MomentumResolSystematicsNumuCC());    

  if (ND::params().GetParameterI("numuCCAnalysis.Systematics.EnableMomScaleOptimization"))
    syst().ReplaceVariationSystematic(baseAnalysis::kMomScale, new MomentumScaleSystematicsNumuCC());

}

//********************************************************************
void numuCCAnalysis::DefineMicroTrees(bool addBase){
//********************************************************************

  // -------- Add variables to the analysis tree ----------------------

  // Variables from baseAnalysis (run, event, ...)
  if (addBase) baseAnalysis::DefineMicroTrees(addBase);


  // --- muon candidate truth variables -------
  AddVarF(  output(), selmu_truemom,    "muon candidate true momentum");
  AddVar4VF(output(), selmu_truepos,    "muon candidate true position");         
  AddVar4VF(output(), selmu_trueendpos, "muon candidate true end position");      
  AddVar3VF(output(), selmu_truedir,    "muon candidate true direction");         

  // --- true muon truth variables -------
  AddVarF(output(),   truemu_truemom,     "true muon true momentum");
  AddVarF(output(),   truemu_truecostheta,"true muon true cos(theta)");

  // --- muon candidate recon variables -----
  AddVar3VF(output(), selmu_dir,        "muon candidate reconstructed direction");
  AddVar4VF(output(), selmu_pos,        "muon candidate reconstructed position");
  AddVar4VF(output(), selmu_endpos,     "muon candidate reconstructed end position");
  AddVarF(  output(), selmu_costheta,   "muon candidate reconstructed cos(theta)");
  AddVarI(  output(), selmu_closest_tpc,"muon candidate closest TPC");
  AddVarI(  output(), selmu_detectors,  "muon candidate detectors");

  AddToyVarF(output(), selmu_mom,      "muon candidate reconstructed momentum");
  AddToyVarF(output(), selmu_theta,    "muon candidate reconstructed polar angle");

  if (ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
    AddToyVarF(output(), selmu_amom,     "muon candidate alternate reconstructed momentum ");
    AddToyVarF(output(), selmu_likemu,   "muon candidate muon TPC PID likelihood"); 
    AddToyVarF(output(), selmu_likemip,  "muon candidate MIP TPC PID likelihood"); 
    AddToyVarF(output(), shmn_mom,       "second highest momentum negative track reconstructed momentum ");  
  }
  else{
    AddVarF(output(), selmu_amom,     "muon candidate alternate reconstructed momentum ");
    AddVarF(output(), selmu_likemu,   "muon candidate muon TPC PID likelihood"); 
    AddVarF(output(), selmu_likemip,  "muon candidate MIP TPC PID likelihood"); 
    AddVarF(output(), shmn_mom,       "second highest momentum negative track reconstructed momentum ");  
  }

  // --- info by tpc
  AddVarVI(output(), selmu_tpc_det,         "muon candidate TPC number",                               selmu_ntpcs);
  AddVarVI(output(), selmu_tpc_nhits,       "muon candidate #hits in each TPC",                        selmu_ntpcs);
  AddVarVI(output(), selmu_tpc_nnodes,      "muon candidate #nodes in each TPC",                       selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_charge,      "muon candidate reconstructed charge in each TPC",         selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_mom,         "muon candidate reconstructed momentum in each TPC",       selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_bfield_mom,  "muon candidate reconstructed momentum with BField refit in each TPC",       selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_efield_mom,  "muon candidate reconstructed momentum with EFiled refit in each TPC",       selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_emom,        "muon candidate reconstructed momentum error in each TPC", selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_truemom,     "muon candidate true momentum in each TPC",                selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_dedx_raw,    "muon candidate raw dEdx (CT) in each TPC",                selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_dedx_expmu,  "muon candidate expected dEdx for muons in each TPC",      selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_dedx_exppi,  "muon candidate expected dEdx for pions in each TPC",      selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_dedx_expele, "muon candidate expected dEdx for electrons in each TPC",  selmu_ntpcs);
  AddVarVF(output(), selmu_tpc_dedx_expp,   "muon candidate expected dEdx for protons in each TPC",    selmu_ntpcs);

  if (ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
    AddToyVarVF(output(), selmu_tpc_pullmu,  "muon candidate muon pull in each TPC",     NMAXTPCS);
    AddToyVarVF(output(), selmu_tpc_pullpi,  "muon candidate pion pull in each TPC",     NMAXTPCS);
    AddToyVarVF(output(), selmu_tpc_pullele, "muon candidate electron pull in each TPC", NMAXTPCS);
    AddToyVarVF(output(), selmu_tpc_pullp,   "muon candidate proton pull in each TPC",   NMAXTPCS);
    AddToyVarVF(output(), selmu_tpc_dedx,    "muon candidate dEdx (CT) in each TPC",     NMAXTPCS);
  }
  else{
    AddVarVF(output(), selmu_tpc_pullmu,  "muon candidate muon pull in each TPC",     selmu_ntpcs);
    AddVarVF(output(), selmu_tpc_pullpi,  "muon candidate pion pull in each TPC",     selmu_ntpcs);
    AddVarVF(output(), selmu_tpc_pullele, "muon candidate electron pull in each TPC", selmu_ntpcs);
    AddVarVF(output(), selmu_tpc_pullp,   "muon candidate proton pull in each TPC",   selmu_ntpcs);
    AddVarVF(output(), selmu_tpc_dedx,    "muon candidate dEdx (CT) in each TPC",     selmu_ntpcs);
  }

  // --- info by FGD
  AddVarVI(output(), selmu_fgd_det,   "muon candidate FGD number",                       selmu_nfgds);
  AddVarVF(output(), selmu_fgd_x,     "muon candidate track length in each FGD",         selmu_nfgds);
  AddVarVF(output(), selmu_fgd_V00,   "muon candidate V00 vertex activity in each FGD",  selmu_nfgds);
  AddVarVF(output(), selmu_fgd_V11,   "muon candidate V11 vertex activity in each FGD",  selmu_nfgds);
  AddVarVF(output(), selmu_fgd_V33,   "muon candidate V33 vertex activity in each FGD",  selmu_nfgds);
  AddVarVF(output(), selmu_fgd_V55,   "muon candidate V55 vertex activity in each FGD",  selmu_nfgds);
  AddVarVF(output(), selmu_fgd_V77,   "muon candidate V77 vertex activity in each FGD",  selmu_nfgds);

  if (ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
    AddToyVarVF(output(), selmu_fgd_pullmu,"muon candidate muon pull in each FGD ", NMAXFGDS);
  }
  else{
    AddVarVF(output(), selmu_fgd_pullmu,"muon candidate muon pull in each FGD ", selmu_nfgds);
  }

  // --- info by ECAL
  
  AddVarVI(output(), selmu_ecal_det,     "The number of hits in the selected track's ECal segment", selmu_necals);

  AddVarVI(output(), selmu_ecal_nhits,     "The number of hits in the selected track's ECal segment", selmu_necals);
  AddVarVI(output(), selmu_ecal_nnodes,     "", selmu_necals);
  AddVarVF(output(), selmu_ecal_length,     "", selmu_necals);
  AddVarMF(output(), selmu_ecal_showerstartpos, "",selmu_necals,-30,4);
  AddVarMF(output(), selmu_ecal_showerendpos, "",selmu_necals,-30,4);
  AddVarMF(output(), selmu_ecal_showerstartdir, "",selmu_necals,-30,3);
  AddVarMF(output(), selmu_ecal_showerenddir, "",selmu_necals,-30,3);

  AddVarVF(output(), selmu_ecal_EMenergy,     "", selmu_necals);
  AddVarVF(output(), selmu_ecal_edeposit,     "", selmu_necals);
  AddVarVI(output(), selmu_ecal_IsShower,     "", selmu_necals);

  AddVarVF(output(), selmu_ecal_mipem,     "MipEm value of the selected track's ECal segment. Negative means more MIP-like, positive means more EM shower-like", selmu_necals);
  AddVarVF(output(), selmu_ecal_mippion,   "MipPion value the selected track's ECal segment. Negative means more MIP-like, positive means more showering pion-like", selmu_necals);
  AddVarVF(output(), selmu_ecal_emhip,     "EmHip value of the selected track's ECal segment. Negative means more EM shower-like, positive means more HIP-like.", selmu_necals);
  AddVarVF(output(), selmu_ecal_containment,     "", selmu_necals);
  AddVarVF(output(), selmu_ecal_trshval,     "", selmu_necals);
  AddVarVI(output(), selmu_ecal_mostupstreamlayerhit,     "", selmu_necals);

}

//********************************************************************
void numuCCAnalysis::DefineTruthTree(){
//********************************************************************

  // Variables from baseAnalysis (run, event, ...)
  baseAnalysis::DefineTruthTree();

  //--- muon variables -------
  //  AddVarF(output(), "truth","truemu_truerec",  "I");  // Is the recon muon a true muon
  AddVarF(  output(), truemu_truemom,      "true muon true momentum");
  AddVarF(  output(), truemu_truecostheta, "true muon true cos(theta)");
  AddVar4VF(output(), truemu_truepos,      "true muon true position");
  AddVar3VF(output(), truemu_truedir,      "true muon true direction");

  AddVarI(output(),fgd2reaction,      "fgd2reaction type");
  AddVarI(output(),fgd2topology,      "fgd2topology type");
  AddVarI(output(),fgd2mectopology,   "fgd2mectopology type");

}

//********************************************************************
bool numuCCAnalysis::Process(){
//********************************************************************

  //------ This function is called for each bunch and each toy experiment  ----------

  // FGD1 or FGD2 analysis
  if (_whichFGD==1)
    return sel().GetSelection("kTrackerNumuCC")->Apply(*_event,box());
  else if (_whichFGD==2)
    return sel().GetSelection("kTrackerFGD2CC")->Apply(*_event,box());
//  else if (_whichFGD==3)
//    return sel().GetSelection("kTrackerFGDsCC")->Apply(*_event,box());

  return false;
}

//********************************************************************
void numuCCAnalysis::FillMicroTrees(bool addBase){
//********************************************************************

  // Variables from baseAnalysis (run, event, ...)
  if (addBase) baseAnalysis::FillMicroTreesBase(addBase);

  // Muon candidate variables
  if (box().HMNtrack){
    if(box().HMNtrack->TrueTrack) {
      AnaTrueVertex *vtx = static_cast<AnaTrueVertex*>(box().HMNtrack->TrueTrack->TrueVertex);
      if(vtx) {
        output().FillVar(truemu_truemom,vtx->LeptonMom);
        output().FillVar(truemu_truecostheta,(Float_t)cos(utils::ArrayToTVector3(vtx->LeptonDir).Angle(utils::ArrayToTVector3(vtx->NuDir))));
      }
    }

    output().FillVar(selmu_costheta,              box().HMNtrack->DirectionStart[2]);
    output().FillVectorVarFromArray(selmu_pos,    box().HMNtrack->PositionStart, 4);
    output().FillVectorVarFromArray(selmu_endpos, box().HMNtrack->PositionEnd, 4);
    output().FillVectorVarFromArray(selmu_dir,    box().HMNtrack->DirectionStart, 3);

    if (!ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
      output().FillVar(selmu_amom, (static_cast<AnaTrack*>(box().HMNtrack))->MomentumMuon);
      // PID likelihoods
      output().FillVar(selmu_likemu, anaUtils::GetPIDLikelihood( *(box().HMNtrack),0)); 
      output().FillVar(selmu_likemip,anaUtils::GetPIDLikelihoodMIP( *(box().HMNtrack))); 
    }

    // Properties of the true particle associated to the muon candidate
    if(  box().HMNtrack->TrueTrack ) {
      output().FillVar(selmu_truemom,         box().HMNtrack->TrueTrack->Momentum);
      output().FillVectorVarFromArray(selmu_truepos,   box().HMNtrack->TrueTrack->Position, 4);
      output().FillVectorVarFromArray(selmu_trueendpos,static_cast<AnaTrueTrack*>(box().HMNtrack->TrueTrack)->PositionEnd, 4);
      output().FillVectorVarFromArray(selmu_truedir,   box().HMNtrack->TrueTrack->Direction, 3);       
    }

    // The closest TPC
    SubDetId::SubDetEnum tpc = anaUtils::GetClosestTPC(*(box().HMNtrack));
    output().FillVar(selmu_closest_tpc, tpc+1);      

    // Track composition
    output().FillVar(selmu_detectors, static_cast<AnaTrack*>(box().HMNtrack)->Detectors);      

    // Info in all TPCs 

    for (Int_t subdet = 0; subdet<3; subdet++) {
      if (!SubDetId::GetDetectorUsed(box().HMNtrack->Detector, static_cast<SubDetId::SubDetEnum >(subdet+2))) continue;
      AnaTpcTrack* TPCSegment = static_cast<AnaTpcTrack*>(anaUtils::GetSegmentInDet( *box().HMNtrack,static_cast<SubDetId::SubDetEnum >(subdet+2)));
      if (!TPCSegment) continue;
      output().FillVectorVar(selmu_tpc_det,     subdet);
      output().FillVectorVar(selmu_tpc_mom,     TPCSegment->Momentum);

      output().FillVectorVar(selmu_tpc_bfield_mom, TPCSegment->RefitMomentum);
      output().FillVectorVar(selmu_tpc_efield_mom, TPCSegment->EFieldRefitMomentum);
      output().FillVectorVar(selmu_tpc_charge,  TPCSegment->Charge);
      output().FillVectorVar(selmu_tpc_nhits,   TPCSegment->NHits);
      output().FillVectorVar(selmu_tpc_nnodes,  TPCSegment->NNodes);
      output().FillVectorVar(selmu_tpc_emom,    TPCSegment->MomentumError);

      if (TPCSegment->TrueTrack)
	output().FillVectorVar(selmu_tpc_truemom, TPCSegment->TrueTrack->Momentum);

      if (TPCSegment->Original->Original)
	output().FillVectorVar(selmu_tpc_dedx_raw,    (static_cast< const AnaTpcTrack* >(TPCSegment->Original->Original))->dEdxMeas);

      output().FillVectorVar(selmu_tpc_dedx_expmu,  TPCSegment->dEdxexpMuon);
      output().FillVectorVar(selmu_tpc_dedx_exppi,  TPCSegment->dEdxexpPion);
      output().FillVectorVar(selmu_tpc_dedx_expele, TPCSegment->dEdxexpEle);
      output().FillVectorVar(selmu_tpc_dedx_expp,   TPCSegment->dEdxexpProton);

      if (!ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
	output().FillVectorVar(selmu_tpc_dedx,   TPCSegment->dEdxMeas);
	output().FillVectorVar(selmu_tpc_pullmu, TPCSegment->Pullmu);
	output().FillVectorVar(selmu_tpc_pullpi, TPCSegment->Pullpi);
	output().FillVectorVar(selmu_tpc_pullele,TPCSegment->Pullele);
	output().FillVectorVar(selmu_tpc_pullp,  TPCSegment->Pullp);    
      }

      output().IncrementCounterForVar(selmu_tpc_det);
    }

    // Info in all FGDs 
    for (Int_t subdet = 0; subdet<2; subdet++) {
      if (!SubDetId::GetDetectorUsed(box().HMNtrack->Detector, static_cast<SubDetId::SubDetEnum >(subdet))) continue;
	  AnaFgdTrack* FGDSegment = static_cast<AnaFgdTrack*>(anaUtils::GetSegmentInDet( *box().HMNtrack,static_cast<SubDetId::SubDetEnum >(subdet)));
      if (!FGDSegment) continue;
      
      //							      output().FillVectorVar("selmu_fgd_det",     (Int_t)FGDSegment->Detector);
      output().FillVectorVar(selmu_fgd_det,      subdet);
      output().FillVectorVar(selmu_fgd_x,        FGDSegment->X);
      output().FillVectorVar(selmu_fgd_V00,      FGDSegment->Vertex0x0);
      output().FillVectorVar(selmu_fgd_V11,      FGDSegment->Vertex1x1);
      output().FillVectorVar(selmu_fgd_V33,      FGDSegment->Vertex3x3);
      output().FillVectorVar(selmu_fgd_V55,      FGDSegment->Vertex5x5);
      output().FillVectorVar(selmu_fgd_V77,      FGDSegment->Vertex7x7);

      if (!ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars"))
	output().FillVectorVar(selmu_fgd_pullmu,   FGDSegment->Pullmu);

      output().IncrementCounterForVar(selmu_fgd_det);
    }

    // Info in all ECALs
    for (Int_t subdet = 0; subdet<9; subdet++) {
      if (!SubDetId::GetDetectorUsed(box().HMNtrack->Detector, static_cast<SubDetId::SubDetEnum >(subdet+6))) continue;
	  AnaEcalTrack* ECALSegment = static_cast<AnaEcalTrack*>(anaUtils::GetSegmentInDet( *box().HMNtrack,static_cast<SubDetId::SubDetEnum >(subdet+6)));
      if (!ECALSegment) continue;

      output().FillVectorVar(selmu_ecal_det,      subdet);
      output().FillVectorVar(selmu_ecal_nhits,      ECALSegment->NHits);
      output().FillVectorVar(selmu_ecal_nnodes,      ECALSegment->NNodes);
      output().FillVectorVar(selmu_ecal_length,      ECALSegment->Length);
      output().FillMatrixVarFromArray(selmu_ecal_showerstartpos, ECALSegment->PositionStart, 4 );
      output().FillMatrixVarFromArray(selmu_ecal_showerendpos, ECALSegment->PositionEnd, 4 );
      output().FillMatrixVarFromArray(selmu_ecal_showerstartdir, ECALSegment->DirectionStart, 3 );
      output().FillMatrixVarFromArray(selmu_ecal_showerenddir, ECALSegment->DirectionEnd, 3 );

      output().FillVectorVar(selmu_ecal_EMenergy,      ECALSegment->EMEnergy);
      output().FillVectorVar(selmu_ecal_edeposit,      ECALSegment->EDeposit);
      output().FillVectorVar(selmu_ecal_IsShower,      ECALSegment->IsShowerLike);

      output().FillVectorVar(selmu_ecal_trshval,      ECALSegment->TrShVal);
      output().FillVectorVar(selmu_ecal_emhip,      ECALSegment->PIDEmHip);
      output().FillVectorVar(selmu_ecal_mippion,      ECALSegment->PIDMipPion);
      output().FillVectorVar(selmu_ecal_mipem,      ECALSegment->PIDMipEm);
      output().FillVectorVar(selmu_ecal_containment,      ECALSegment->Containment);
      output().FillVectorVar(selmu_ecal_mostupstreamlayerhit,      ECALSegment->MostUpStreamLayerHit);
  
      output().IncrementCounterForVar(selmu_ecal_det);
	
    }
  
  }
  
}

//********************************************************************
void numuCCAnalysis::FillToyVarsInMicroTrees(bool addBase){
//********************************************************************

  // Fill all tree variables that vary for each virtual analysis (toy experiment)

  // Fill the common variables
  if (addBase) baseAnalysis::FillToyVarsInMicroTreesBase(addBase);

  // variables specific for this analysis
  if (box().HMNtrack){
    output().FillToyVar(selmu_mom, box().HMNtrack->Momentum);    
    output().FillToyVar(selmu_theta,   (Float_t)acos(box().HMNtrack->DirectionStart[2]));       

    if (ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
      output().FillToyVar(selmu_amom, (static_cast<AnaTrack*>(box().HMNtrack))->MomentumMuon);
      // PID likelihoods
      output().FillToyVar(selmu_likemu, anaUtils::GetPIDLikelihood( *(box().HMNtrack),0)); 
      output().FillToyVar(selmu_likemip,anaUtils::GetPIDLikelihoodMIP( *(box().HMNtrack))); 
      
      // Info in all TPCs 
      for (Int_t subdet = 0; subdet<3; subdet++) {
	if (!SubDetId::GetDetectorUsed(box().HMNtrack->Detector, static_cast<SubDetId::SubDetEnum >(subdet+2))) continue;
	AnaTpcTrack* TPCSegment = static_cast<AnaTpcTrack*>(anaUtils::GetSegmentInDet( *box().HMNtrack, static_cast<SubDetId::SubDetEnum >(subdet+2)));
	if (!TPCSegment) continue;
	output().FillToyVectorVar(selmu_tpc_dedx,   TPCSegment->dEdxMeas, subdet);
	output().FillToyVectorVar(selmu_tpc_pullmu, TPCSegment->Pullmu,   subdet);
	output().FillToyVectorVar(selmu_tpc_pullpi, TPCSegment->Pullpi,   subdet);
	output().FillToyVectorVar(selmu_tpc_pullele,TPCSegment->Pullele,  subdet);
	output().FillToyVectorVar(selmu_tpc_pullp,  TPCSegment->Pullp,    subdet);    
      }

      // Info in all FGDs 
      for (Int_t subdet = 0; subdet<2; subdet++) {
	if (!SubDetId::GetDetectorUsed(box().HMNtrack->Detector, static_cast<SubDetId::SubDetEnum >(subdet))) continue;
	AnaFgdTrack* FGDSegment = static_cast<AnaFgdTrack*>(anaUtils::GetSegmentInDet(*(box().HMNtrack), static_cast<SubDetId::SubDetEnum >(subdet)));
	if (!FGDSegment) continue;
	output().FillToyVectorVar(selmu_fgd_pullmu, FGDSegment->Pullmu, subdet);
      }
    }
  }
  if (ND::params().GetParameterI("numuCCAnalysis.MicroTrees.AdditionalToyVars")){
    if (box().SHMNtrack){
      output().FillToyVar(shmn_mom, box().SHMNtrack->Momentum);    
    }
  }
}

//********************************************************************
bool numuCCAnalysis::CheckFillTruthTree(const AnaTrueVertex& vtx){
//********************************************************************

  bool numuCC=vtx.ReacCode>0 && vtx.ReacCode<30 && vtx.NuPDG==14;// && vtx.LeptonPDG==13;
  SubDetId::SubDetEnum _fgdID;
  if (_whichFGD==1) _fgdID = SubDetId::kFGD1;
  if (_whichFGD==2) _fgdID = SubDetId::kFGD2;
//  if (_whichFGD==3) _fgdID = SubDetId::kFGD;
  return (anaUtils::InFiducialVolume(_fgdID, vtx.Position) && numuCC);
}

//********************************************************************
bool numuCCAnalysis::CheckFillRooTracker(const AnaTrueVertex& vtx){
//********************************************************************

  bool numuCC=vtx.ReacCode>0 && vtx.ReacCode<30 && vtx.NuPDG==14;
  return numuCC;
}

//********************************************************************
void numuCCAnalysis::FillTruthTree(const AnaTrueVertex& vtx){
//********************************************************************

  // Fill the common variables
  baseAnalysis::FillTruthTreeBase(vtx);

  // Fill truth categories for FGD2, same categories with prefix "fgd2", i,e, "fgd2reaction" etc.
  // We could directly set which FGD in the usual categories,
  // but better having dedicated FGD2 categories for the selection with both FGDs
  anaUtils::FillTruthTreeCategories(vtx, "fgd2", SubDetId::kFGD2);
  // to do: automatize the truth categories filling in Analysis loop,
  // as it is for categories in the default tree, so to avoid this specific code for FGD2)
  output().FillVar(fgd2reaction,      ND::categ().GetCode("fgd2reaction"));
  output().FillVar(fgd2topology,      ND::categ().GetCode("fgd2topology"));
  output().FillVar(fgd2mectopology,   ND::categ().GetCode("fgd2mectopology"));

  // Muon true variables
  output().FillVar(truemu_truecostheta,(Float_t)cos(utils::ArrayToTVector3(vtx.LeptonDir).Angle(utils::ArrayToTVector3(vtx.NuDir))));
  output().FillVar(truemu_truemom,vtx.LeptonMom);
  output().FillVectorVarFromArray(truemu_truedir, vtx.LeptonDir,3);
  output().FillVectorVarFromArray(truemu_truepos, vtx.Position, 4);

  // Muon true-rec association.
  //  if (box().HMNtrack)
  //    if (box().HMNtrack->TrueTrack)
  //      output().FillVar("mu_true_rec", (box().HMNtrack->TrueTrack->PDG==13));   // Is the recon muon a true muon ?
}

//********************************************************************
void numuCCAnalysis::FillCategories(){
//********************************************************************

  // Fill the track categories for color drawing

  // For the muon candidate
  anaUtils::FillCategories(_event, static_cast<AnaTrack*>(box().HMNtrack),"");


  // For FGD2, same categories with prefix "fgd2", i,e, "fgd2reaction" etc.
  // We could directly set which FGD in the usual categories,
  // but better having dedicated FGD2 categories for the selection with both FGDs
  anaUtils::FillCategories(_event, static_cast<AnaTrack*>(box().HMNtrack), "fgd2", SubDetId::kFGD2);
}
